      Subroutine Fssi_Mc_Jifen2(Stress,Statev,Ddsdde,Sse,Spd,Scd,Rpl,Ddsddt,Drplde,Drpldt,   &
                                Stran,Dstran,Time,Dtime,Temp,Dtemp,Predef,Dpred,Cmname,   &
                                Ndi,Nshr,Ntens,Nstatv,Props,Nprops,Coords,Drot,Pnewdt,   &
                                Celent,Dfgrd0,Dfgrd1,Noel,Npt,Layer,Kspt,Jstep,Kinc,Func)
!
!     Include 'Aba_Param.Inc'    ! Linux不支持Include
      Implicit Real*8(A-H,O-Z)
      Parameter (Nprecd=2)
!
      Character*80 Cmname
      Dimension Stress(Ntens),Statev(Nstatv),Ddsdde(Ntens,Ntens),Ddsddt(Ntens),Drplde(Ntens),   &
                Stran(Ntens),Dstran(Ntens),Time(2),Predef(1),Dpred(1),Props(Nprops),Coords(3),   &
                Drot(3,3),Dfgrd0(3,3),Dfgrd1(3,3),Jstep(4)
      Integer::Func,Ntest
      Dimension Ps(3)
      Parameter (Zeropointfive=0.5D0,One=1.0D0,Two=2.0D0,Three=3.0D0,Six=6.0D0)
      Character*5      :: Name


      Double Precision::Param(Nstatv),Dstress(Ntens)

!*****************************************************
!    ** 对应力和应变进行换号，本构中 应力以压为正
!                                    应变以压缩为正
      Stress=Stress
      Dstran=Dstran

!*****************************************************

      Dstress=0.0D0
      Ntest=4
      
      If (Int(Func).Eq.Int(1))Then
         Call Getstress(Props,Nprops,Statev,Nstatv,Stran,Dstran,Stress,Dstress,Ntens,Ddsdde,Noel,Npt,Nstatv+Ntens,Jstep,Ntest)
      Else
         Call Getdmatrx(Props,Nprops,Statev,Nstatv,Stran,Dstran,Stress,Dstress,Ddsdde,Ntens,Noel,Npt,Jstep,Ntest)
      End If

!*****************************************************
!    ** 对应力和应变进行换号，本构中 应力以压为正
!                                    应变以压缩为正
      Stress=Stress
      Dstran=Dstran

!*****************************************************

      Return
      End


!********************************************************************************
!  One          *************** Get Stress *************************
!********************************************************************************
      Subroutine Getstress(Propd,Lprpd,Param,Lpara,Strain,Dstan,Stress,Dstre,   &
                           Nstre,Dmatx,Ielem,Igaus,Mgmod,Ktest,Ntest)

!     Subroutine To Integrate The Incremental Stresses
!.....增量应力的积分
!....
      Implicit None
!.... Input Parameters:
      Integer Ielem,Igaus,Model,Lprpd,Iswdp,Icdat,Icout,Mgmod
      Integer Nstre,Lpara,Mdivt,Icprt,Ichpc,Ntest
      Double Precision Propd(Lprpd),Voidr,Dstan(Nstre),Tolst,Tbeta
!.... Input/Output Parameters:
      Double Precision Gausm(Mgmod),Porep,Dpore,Gpcod(3),Cord0(3),Almin
!.... Output Parameters:
      Double Precision Dstre(Nstre)
!.... Passing Variables:
      Double Precision Dmatx(Nstre,Nstre)
      Double Precision Tolcl1_Umat,Tolcl2_Umat,Tolcl3_Umat
!.... Local Variables:
      Double Precision Stora(Mgmod),Storb(Mgmod),Strea(Nstre)
      Double Precision Storc(Mgmod),Strec(Nstre)
      Double Precision Streb(Nstre),Avalu(6)
      Double Precision Stana(Nstre),Cons1,Cons2,Total,Alpha,Tempa
      Integer Ndivt,Idivt,Istre,Ktest(Ntest),Kdivt,Itest,Ksubv
!.... Constants
      Double Precision Rep99,Real0,Real1
      Integer Kprin
      Logical Debua,Debug
      Parameter (Debua=.True.,Debug=.False.)
      Logical D0Eq
      External D0Eq
!.... Printing Keys
      Integer Kprt0,Kprt1,Kprt2,Kprt3,Kprt4
      Integer Kprt5,Kprt6,Krpt7,Krpt8,Kprt9
!.... Initialization

!     Local Variable
      Integer::I
      Double Precision:: Stress(Nstre),Param(Lpara),Tolcl1,Strain(Nstre)
      Double Precision:: Ndepv

      Double Precision:: Alpha0,Alpha1,Almax
      Double Precision:: Temp0(Mgmod),Temp1(Mgmod),Temp2(Mgmod)
      Double Precision:: Temp3(Mgmod),Temp4(Mgmod)
      Double Precision:: Temp1_Dsh(Mgmod),Temp2_Dsh(Mgmod)
      Double Precision:: Temp3_Dsh(Mgmod),Temp4_Dsh(Mgmod)
      Double Precision:: Dstre1(Nstre),Dstre2(Nstre)
      Double Precision:: Dstre3(Nstre),Dstre4(Nstre)
      Double Precision:: Temprk2(Mgmod),Temprk3(Mgmod),Temprk4(Mgmod)
      Double Precision::Dstre5(Mgmod),Dstre6(Mgmod)
      Double Precision::Temp5(Mgmod),Temp6(Mgmod)
      Double Precision::Temp5_Dsh(Mgmod),Temp6_Dsh(Mgmod)
      Double Precision::Temprk6_Bar(Mgmod),Temprk6(Mgmod)


      Rep99=0.99D0
      Real0=0.0D0
      Real1=1.0D0
      Kprin=100

      Tbeta=0.95D0                   !读入应变细分的Beta值         !Beta Value For Automatic Strain Subdivision

      Model=1
      Total=Real0
      Alpha=Real1
      Cons2=Tolst*Tbeta
      Cons1=Cons2
      Iswdp=2
      Cord0=0.0D0
      Avalu=0.0D0

      Dstre=0.0D0

!     Mgmod=Lpara+Nstre

      Do I=1,Nstre
         Gausm(I)=Stress(I)
      End Do

      Do I=1,Lpara
         Gausm(Nstre+I)=Param(I)
      End Do

      Kdivt=0

      Ksubv=Propd(7)              !选择积分算法类型取值1-7
      Tolst=Propd(8)              !读入非线性收敛的精度Tolst    !Tolerance For Strain Subdivision:
      Mdivt=Propd(9)              !读入应变细分的最大数量Mdivt  !Maximum Number Of Strain Subdivision  
      Almin=Propd(10)             !读入应变细分的最小步长Almin  !Minimum Size Of Strain Subdivision:
      
      Goto (1100,1200,1300,1400,1500),Ksubv
1100  Continue            ! 不细分
      Ktest(1)=1
      Ndepv=Ndepv+1
      Call Dcopy_Umat (Mgmod,Gausm,1,Stora,1)
      Call Consti_Umat(Model,Propd,Lprpd,Iswdp,Icdat,Icout,Dstre,Dstan,Dmatx,Gausm,Mgmod,Ielem,   &
                       Igaus,Voidr,Nstre,Lpara,Icprt,Ktest,Ntest,Porep,Dpore,Gpcod,Ichpc,Avalu,Cord0)

      Goto 1900

1200  Continue            ! 给定最大细分数（>1）
      Ktest(1)=1
!.... Store Away The Initial State
      Call Dcopy_Umat (Mgmod,Gausm,1,Stora,1)

!.... Calculate The Stress Increment With No Subdvision
      Ndepv=Ndepv+1
      Call Consti_Umat(Model,Propd,Lprpd,Iswdp,Icdat,Icout,Dstre,Dstan,Dmatx,Gausm,Mgmod,Ielem,   &
                       Igaus,Voidr,Nstre,Lpara,Icprt,Ktest,Ntest,Porep,Dpore,Gpcod,Ichpc,Avalu,Cord0)
!.... Calculate The Stress Increment Based On The Final State
      Ktest(1)=0
      Call Dcopy_Umat (Mgmod,Gausm,1,Storb,1)

      Ndepv=Ndepv+1
      Call Consti_Umat(Model,Propd,Lprpd,Iswdp,Icdat,Icout,Streb,Dstan,Dmatx,Storb,Mgmod,Ielem,   &
                       Igaus,Voidr,Nstre,Lpara,Icprt,Ktest,Ntest,Porep,Dpore,Gpcod,Ichpc,Avalu,Cord0)
!.... Compare The Calculate The Number Of Subdivision
      Cons1=Tolcl1_Umat(Dstre,Streb,Nstre,Stora,Tolst)
      Ndivt=Cons1/Tolst+Rep99
      If (Ndivt.Gt.Mdivt) Ktest(3)=Ktest(3)+1

      Ndivt=Min(Max(Ndivt,1),Mdivt)
!.... Subdivision Is Needed
      If (Ndivt.Gt.1) Then
!.... Calculate The Subdividing Strain
         Ktest(1)=2
         Do Istre=1,Nstre
            Stana(Istre)=Dstan(Istre)/Ndivt
            Dstre(Istre)=Real0
            Strea(Istre)=Real0
         End Do
!.... Restore The Initial State
         Call Dcopy_Umat(Mgmod,Stora,1,Gausm,1)

         Do Idivt=1,Ndivt
            Ndepv=Ndepv+1
            Call Consti_Umat(Model,Propd,Lprpd,Iswdp,Icdat,Icout,Strea,Stana,Dmatx,Gausm,Mgmod,Ielem,   &
                             Igaus,Voidr,Nstre,Lpara,Icprt,Ktest,Ntest,Porep,Dpore,Gpcod,Ichpc,Avalu,Cord0)
            Do Istre=1,Nstre
               Dstre(Istre)=Dstre(Istre)+Strea(Istre)
            End Do
            Strea=Real0
         End Do
      Else
         Dstre=(Streb+Dstre)*0.5D0
      End If
      Goto 1900

1300  Continue            ! 给定最小增量应变（>almin）
      Ktest(1)=1
!.... Store Away The Initial State
      Call Dcopy_Umat (Mgmod,Gausm,1,Stora,1)

!.... Calculate The Stress Increment With No Subdvision
      Ndepv=Ndepv+1
      Call Consti_Umat(Model,Propd,Lprpd,Iswdp,Icdat,Icout,Dstre,Dstan,Dmatx,Gausm,Mgmod,Ielem,   &
                       Igaus,Voidr,Nstre,Lpara,Icprt,Ktest,Ntest,Porep,Dpore,Gpcod,Ichpc,Avalu,Cord0)
!.... Calculate The Stress Increment Based On The Final State
      Ktest(1)=0
      Call Dcopy_Umat (Mgmod,Gausm,1,Storb,1)

      Ndepv=Ndepv+1
      Call Consti_Umat(Model,Propd,Lprpd,Iswdp,Icdat,Icout,Streb,Dstan,Dmatx,Storb,Mgmod,Ielem,   &
                       Igaus,Voidr,Nstre,Lpara,Icprt,Ktest,Ntest,Porep,Dpore,Gpcod,Ichpc,Avalu,Cord0)
!.... Compare The Calculate The Number Of Subdivision
      Cons1=Tolcl2_Umat(Dstre,Streb,Nstre,Stora,Tolst)

      Ndivt=1
      Ndivt=Max(Max(Ndivt,1),Int(0.05*Mdivt))
!.... Subdivision Is Needed
      If (Cons1.Gt.Tolst) Then
!.... Calculate The Subdividing Strain
         Ktest(1)=2
         Do Istre=1,Nstre
            Stana(Istre)=Dstan(Istre)/Ndivt
            Dstre(Istre)=Real0
            Strea(Istre)=Real0
         End Do
!.... Restore The Initial State
         Call Dcopy_Umat(Mgmod,Stora,1,Gausm,1)

         Do Idivt=1,Ndivt
            Ndepv=Ndepv+1
            Call Consti_Umat(Model,Propd,Lprpd,Iswdp,Icdat,Icout,Strea,Stana,Dmatx,Gausm,Mgmod,Ielem,   &
                             Igaus,Voidr,Nstre,Lpara,Icprt,Ktest,Ntest,Porep,Dpore,Gpcod,Ichpc,Avalu,Cord0)
            Do Istre=1,Nstre
               Dstre(Istre)=Dstre(Istre)+Strea(Istre)
            End Do
            Strea=Real0
         End Do
      Else
         Dstre=(Streb+Dstre)*0.5D0
      End If
      Goto 1900

!************************************************************************************
!****************                   三阶龙格库塔                *********************
!************************************************************************************
1400  Continue
      Ktest(1)=1
!.... Store Away The Initial State
      Call Dcopy_Umat(Mgmod,Gausm,1,Stora,1)

!.... Calculate The Stress Increment With No Subdvision
      Ndepv=Ndepv+1
      Call Consti_Umat(Model,Propd,Lprpd,Iswdp,Icdat,Icout,Dstre,Dstan,Dmatx,Gausm,Mgmod,Ielem,   &
                       Igaus,Voidr,Nstre,Lpara,Icprt,Ktest,Ntest,Porep,Dpore,Gpcod,Ichpc,Avalu,Cord0)
!.... Calculate The Stress Increment Based On The Final State
      Ktest(1)=0
      Call Dcopy_Umat(Mgmod,Gausm,1,Storb,1)

      Ndepv=Ndepv+1
      Call Consti_Umat(Model,Propd,Lprpd,Iswdp,Icdat,Icout,Streb,Dstan,Dmatx,Storb,Mgmod,Ielem,   &
                       Igaus,Voidr,Nstre,Lpara,Icprt,Ktest,Ntest,Porep,Dpore,Gpcod,Ichpc,Avalu,Cord0)
!.... Compare The Calculate The Number Of Subdivision
      Cons1=Tolcl3_Umat(Gausm,Storb,Nstre,Lpara,Stora,Tolst,1)
!....    如果前后两步的应力增量值完全一样，则认为可以是弹性。但是这里面还有一个新的问题，
!        如果是理想弹塑性，那么每一步的值都一样，因此，前后两步骤之间的差值也是0
!        为解决这问题，将硬化模量以及其他塑性应变的值也纳入考虑，因为在进入理想弹塑性的时候，
!        前后两步的塑性应变绝对不一样，所以将传入，Tolcl3_Umat的量换成应力和硬化参量的集合
!.... Subdivision Is Needed
      If (Cons1.Le.Tolst) Then
!说明当前应力增量步全是弹性
         Dstre=(Streb+Dstre)*0.5D0
!           Go To 1430
      Else
!说明存在塑性
!Step 1：初始化

         Total=Real0
         Alpha0=Real1            !步长控制因子
         Alpha1=Real1
         Cons2=Tolst*Tbeta
         Cons1=Cons2
         Dstre=Real0
         Kdivt=0                 !记录计算次数
         Almax=2.0D0

         Temp0=Real0
         Temp1=Real0
         Temp2=Real0
         Temp3=Real0
         Temp4=Real0
         Temprk2=Real0
         Temprk3=Real0

         Temp1_Dsh=Real0
         Temp2_Dsh=Real0
         Temp3_Dsh=Real0
         Temp4_Dsh=Real0

         Dstre1(Nstre)=Real0
         Dstre2(Nstre)=Real0
         Dstre3(Nstre)=Real0
         Dstre4(Nstre)=Real0

!保存最初始的状态结果
1401     Call Dcopy_Umat(Mgmod,Stora,1,Gausm,1)

1402     Call Dcopy_Umat(Mgmod,Gausm,1,Temp0,1)     !将上一个增量步的计算结果保存到 Temp0 中

1410     Continue
!Step 2：进行应变分割，确定子步的步长 Alpha

         Alpha1=Max(Alpha0,Almin) !Almin 是最小的时间步长系数
         Alpha =Min(Alpha1,Almax) !Almax 是最大的时间步长系数,允许大两倍
         Tempa=Real1-Total        !增量时间步长

!判断最后增量步的大小。其根本目的在于要走完全部的应变增量
         If(Alpha.Ge.Tempa) Then
            Alpha=Tempa
         End If

         Call Dcopy_Umat(Mgmod,Temp0,1,Temp1,1)                !将上一个增量步的计算结果保存到 Temp1 中

         Do Istre=1,Nstre
            Stana(Istre)=Alpha*Dstan(Istre)
         End Do
!Step3_1: 计算一阶增量
         Call Consti_Umat(Model,Propd,Lprpd,Iswdp,Icdat,Icout,Dstre1,Stana,Dmatx,Temp1,Mgmod,Ielem,   &
                          Igaus,Voidr,Nstre,Lpara,Icprt,Ktest,Ntest,Porep,Dpore,Gpcod,Ichpc,Avalu,Cord0)
         Temp1_Dsh=Temp1-Temp0                              !增量Ds1, Dh1

!Step3_2: 计算二阶增量                            !采用了中点公式
         Temp2=Temp0+Temp1_Dsh*0.5D0
         Call Consti_Umat(Model,Propd,Lprpd,Iswdp,Icdat,Icout,Dstre2,Stana,Dmatx,Temp2,Mgmod,Ielem,   &
                          Igaus,Voidr,Nstre,Lpara,Icprt,Ktest,Ntest,Porep,Dpore,Gpcod,Ichpc,Avalu,Cord0)
         Temp2_Dsh=Temp2-(Temp0+Temp1_Dsh*0.5D0)            !增量Ds2, Dh2

!Step3_3: 计算三阶增量
         Temp3=Temp0-Temp1_Dsh+Temp2_Dsh*2.0D0
         Call Consti_Umat(Model,Propd,Lprpd,Iswdp,Icdat,Icout,Dstre3,Stana,Dmatx,Temp3,Mgmod,Ielem,   &
                          Igaus,Voidr,Nstre,Lpara,Icprt,Ktest,Ntest,Porep,Dpore,Gpcod,Ichpc,Avalu,Cord0)
         Temp3_Dsh=Temp3-(Temp0-Temp1_Dsh+Temp2_Dsh*2.0D0)  !增量Ds3, Dh4

!Step3：计算三阶龙格库塔的计算结果
         Temprk2=0.0D0
         Temprk3=0.0D0
         Temprk2=Temp0+Temp2_Dsh   ! 采用了中点公式
         Temprk3=Temp0+(1.0D0*Temp1_Dsh+4.0D0*Temp2_Dsh+1.0D0*Temp3_Dsh)/6.0D0

!Step4： 计算误差
         Cons2=0.0D0
         Cons2=Tolcl3_Umat(Temprk2,Temprk3,Nstre,Lpara,Temp0,Tolst,2)
         If (Cons2.Lt.Tolst) Then
            Kdivt=Kdivt+1
            If(Kdivt.Gt.Mdivt)Print*, 'Substeps Exceeds: Kdivt=',Kdivt
            Total=Total+Alpha
            If(Total.Eq.1.0D0)Then
                Gausm=Temprk3
                Do Istre=1,Nstre
                   Dstre(Istre)=Gausm(Istre)-Stora(Istre)
                End Do
                Go To 1430
            Else
               Gausm=Temprk3
               Alpha0=1.0005D0*Alpha
               Go To 1402
            End If
         Else
            Kdivt=Kdivt+1
            If(Kdivt.Gt.Mdivt)Print*, 'Substeps Exceeds: Kdivt=',Kdivt
            Alpha0=0.25D0*Alpha
            Go To 1402
         End If
      End If
1430  Continue

      Goto 1900
      
!************************************************************************************
!****************                   四阶龙格库塔                *********************
!************************************************************************************
1500  Continue
      Ktest(1)=1
!.... Store Away The Initial State
      Call Dcopy_Umat(Mgmod,Gausm,1,Stora,1)

!.... Calculate The Stress Increment With No Subdvision
      Ndepv=Ndepv+1
      Call Consti_Umat(Model,Propd,Lprpd,Iswdp,Icdat,Icout,Dstre,Dstan,Dmatx,Gausm,Mgmod,Ielem,   &
                       Igaus,Voidr,Nstre,Lpara,Icprt,Ktest,Ntest,Porep,Dpore,Gpcod,Ichpc,Avalu,Cord0)
!.... Calculate The Stress Increment Based On The Final State
      Ktest(1)=0
      Call Dcopy_Umat(Mgmod,Gausm,1,Storb,1)

      Ndepv=Ndepv+1
      Call Consti_Umat(Model,Propd,Lprpd,Iswdp,Icdat,Icout,Streb,Dstan,Dmatx,Storb,Mgmod,Ielem,   &
                       Igaus,Voidr,Nstre,Lpara,Icprt,Ktest,Ntest,Porep,Dpore,Gpcod,Ichpc,Avalu,Cord0)
!.... Compare The Calculate The Number Of Subdivision
      Cons1=Tolcl3_Umat(Gausm,Storb,Nstre,Lpara,Stora,Tolst,1)
!....    如果前后两步的应力增量值完全一样，则认为可以是弹性。但是这里面还有一个新的问题，
!        如果是理想弹塑性，那么每一步的值都一样，因此，前后两步骤之间的差值也是0
!        为解决这问题，将硬化模量以及其他塑性应变的值也纳入考虑，因为在进入理想弹塑性的时候，
!        前后两步的塑性应变绝对不一样，所以将传入，Tolcl3_Umat的量换成应力和硬化参量的集合
!.... Subdivision Is Needed
      If (Cons1.Le.Tolst) Then
!说明当前应力增量步全是弹性
         Dstre=(Streb+Dstre)*0.5D0
!           Go To 1530
      Else
!说明存在塑性
!Step 1：初始化

         Total=Real0
         Alpha0=Real1            !步长控制因子
         Alpha1=Real1
         Cons2=Tolst*Tbeta
         Cons1=Cons2
         Dstre=Real0
         Kdivt=0                 !记录计算次数
         Almax=2.0D0

         Temp0=Real0
         Temp1=Real0
         Temp2=Real0
         Temp3=Real0
         Temp4=Real0
         Temprk3=Real0
         Temprk4=Real0

         Temp1_Dsh=Real0
         Temp2_Dsh=Real0
         Temp3_Dsh=Real0
         Temp4_Dsh=Real0

         Dstre1=Real0
         Dstre2=Real0
         Dstre3=Real0
         Dstre4=Real0

!保存最初始的状态结果
1501     Call Dcopy_Umat(Mgmod,Stora,1,Gausm,1)

1502     Call Dcopy_Umat(Mgmod,Gausm,1,Temp0,1)     !将上一个增量步的计算结果保存到 Temp0 中

1510     Continue
!Step 2：进行应变分割，确定子步的步长 Alpha

         Alpha1=Max(Alpha0,Almin) !Almin 是最小的时间步长系数
         Alpha =Min(Alpha1,Almax) !Almax 是最大的时间步长系数,允许大两倍
         Tempa=Real1-Total        !增量时间步长

!判断最后增量步的大小。其根本目的在于要走完全部的应变增量
         If(Alpha.Ge.Tempa) Then
            Alpha=Tempa
         End If

!Step3: 计算三阶龙格库塔的结果
         Call Dcopy_Umat(Mgmod,Temp0,1,Temp1,1)                !将上一个增量步的计算结果保存到 Temp1 中

         Do Istre=1,Nstre
            Stana(Istre)=Alpha*Dstan(Istre)
         End Do

!Step3_1: 计算一阶增量
         Call Consti_Umat(Model,Propd,Lprpd,Iswdp,Icdat,Icout,Dstre1,Stana,Dmatx,Temp1,Mgmod,Ielem,   &
                          Igaus,Voidr,Nstre,Lpara,Icprt,Ktest,Ntest,Porep,Dpore,Gpcod,Ichpc,Avalu,Cord0)
         Temp1_Dsh=Temp1-Temp0                              !增量Ds1, Dh1

!Step3_2: 计算二阶增量                            !采用了中点公式
         Temp2=Temp0+Temp1_Dsh*0.5D0
         Call Consti_Umat(Model,Propd,Lprpd,Iswdp,Icdat,Icout,Dstre2,Stana,Dmatx,Temp2,Mgmod,Ielem,   &
                          Igaus,Voidr,Nstre,Lpara,Icprt,Ktest,Ntest,Porep,Dpore,Gpcod,Ichpc,Avalu,Cord0)
         Temp2_Dsh=Temp2-(Temp0+Temp1_Dsh*0.5D0)            !增量Ds2, Dh2

!Step3_3: 计算三阶增量
         Temp3=Temp0-Temp1_Dsh+Temp2_Dsh*2.0D0
         Call Consti_Umat(Model,Propd,Lprpd,Iswdp,Icdat,Icout,Dstre3,Stana,Dmatx,Temp3,Mgmod,Ielem,   &
                          Igaus,Voidr,Nstre,Lpara,Icprt,Ktest,Ntest,Porep,Dpore,Gpcod,Ichpc,Avalu,Cord0)
         Temp3_Dsh=Temp3-(Temp0-Temp1_Dsh+Temp2_Dsh*2.0D0)  !增量Ds3, Dh4

!Step3_4：计算三阶龙格库塔的计算结果
         Temprk3=Temp0+(1.0D0*Temp1_Dsh+4.0D0*Temp2_Dsh+1.0D0*Temp3_Dsh)/6.0D0

!****************************************************************************************************
!Step4: 计算四阶龙格库塔的结果
!因为四阶龙格库塔表达式中，一阶和二阶增量的表达式完全一致
         Temp3=0.0D0
!Step4_3: 计算三阶增量
         Temp3=Temp0+Temp2_Dsh*0.5D0
         Call Consti_Umat(Model,Propd,Lprpd,Iswdp,Icdat,Icout,Dstre3,Stana,Dmatx,Temp3,Mgmod,Ielem,   &
                          Igaus,Voidr,Nstre,Lpara,Icprt,Ktest,Ntest,Porep,Dpore,Gpcod,Ichpc,Avalu,Cord0)
         Temp3_Dsh=Temp3-(Temp0+Temp2_Dsh*0.5D0)  !增量Ds3, Dh4

!Step4_4：计算三阶龙格库塔的计算结果
!Step4_3: 计算三阶增量
         Temp4=Temp0+Temp3_Dsh
         Call Consti_Umat(Model,Propd,Lprpd,Iswdp,Icdat,Icout,Dstre4,Stana,Dmatx,Temp4,Mgmod,Ielem,   &
                          Igaus,Voidr,Nstre,Lpara,Icprt,Ktest,Ntest,Porep,Dpore,Gpcod,Ichpc,Avalu,Cord0)
         Temp4_Dsh=Temp4-(Temp0+Temp3_Dsh)  !增量Ds3, Dh4

!Step4_4：计算三阶龙格库塔的计算结果
         Temprk4=Temp0+(1.0D0*Temp1_Dsh+2.0D0*Temp2_Dsh+2.0D0*Temp3_Dsh+1.0D0*Temp4_Dsh)/6.0D0

!Step4： 计算误差
         Cons2=0.0D0
         Cons2=Tolcl3_Umat(Temprk3,Temprk4,Nstre,Lpara,Temp0,Tolst,2)
         If (Cons2.Lt.Tolst) Then
            Kdivt=Kdivt+1
            If(Kdivt.Gt.Mdivt)Print*, 'Substeps Exceeds: Kdivt=',Kdivt
            Total=Total+Alpha
            If(Total.Eq.1.0D0)Then
               Gausm=Temprk4
               Do Istre=1,Nstre
                  Dstre(Istre)=Gausm(Istre)-Stora(Istre)
               End Do
               Go To 1530
            Else
               Gausm=Temprk4
               Alpha0=1.0005D0*Alpha
               Go To 1502
            Endif
         Else
            Kdivt=Kdivt+1
            If(Kdivt.Gt.Mdivt)Print*, 'Substeps Exceeds: Kdivt=',Kdivt
            Alpha0=0.25D0*Alpha
            Go To 1502
         End If
      End If
1530  Continue

      Goto 1900

1900  Continue

      Do I=1,Nstre
         Stress(I)=Stora(I)+Dstre(I)
      End Do
      Do I=1,Lpara
         Param(I)=Gausm(Nstre+I)
      End Do


      Return
      End Subroutine

!*************************************************************************
!*************   得到本构中的应力解 *******

      Subroutine Consti_Umat(Model,Propd,Lprpd,Iswdp,Icdat,Icout,Dstre,Dstan,Dmatx,Gausm,Mgmod,Ielem,   &
                             Igaus,Voidr,Nstre,Lpara,Icprt,Ktest,Ntest,Porep,Dpore,Gpcod,Ichpc,Avalu,Cord0)
      Implicit None
      Integer Model,Lprpd,Iswdp,Icdat,Icout,Mgmod,Ielem,Igaus,Nstre
      Double Precision Porep,Dpore,Gpcod(3),Cord0(3)
      Integer Lpara,Icprt,Ichpc
      Integer Ntest,Ktest(Ntest)
      Double Precision Propd(Lprpd),Voidr,Avalu(6)
      Double Precision Dstre(Nstre),Dstan(Nstre),Dmatx(Nstre,Nstre)
      Double Precision Gausm(Mgmod)
      Integer Mstre,Istre
      Parameter (Mstre=6)
      Double Precision Dstanx(Mstre),Voidx
      Double Precision Fstre(Mstre)
      Logical Super
      Integer Isupr,Itemp,Ischn,I
      Parameter (Super=.False.,Ischn=89)
      Double Precision Dtemp
      Double Precision Vsmal
      Parameter (Vsmal=1.0D-30)

!     Local Variable
      Double Precision Stress(Nstre),Param(Lpara)


      Do I=1,Nstre
         Stress(I)=Gausm(I)
      End Do
      Do I=1,Lpara
         Param(I)=Gausm(Nstre+I)
      End Do

      Call Mcoul3_Umat(Propd,Lprpd,Iswdp,Icdat,Icout,Dstre,Dstan,Dmatx,Stress,   &
                       Param,Lpara,Ielem,Igaus,Voidr,Ktest,Ntest,Nstre)
      Voidr=Voidx
      
      Do I=1,Nstre
         Gausm(I)=Stress(I)
      End Do
      Do I=1,Lpara
         Gausm(Nstre+I)=Param(I)
      End Do

      Return
      End Subroutine


!*************************************************************************
!          *************   Function And Subroutine *******

      Function Tolcl1_Umat(Strea,Streb,Nstre,Stora,Tolst)
!....
!     Incremental Stress Norm Function 1
!....
!.... Input Parameters:
      Integer Nstre
      Double Precision Strea(Nstre),Streb(Nstre)
      Double Precision Tolst,Stora(Nstre)

!.... Output Parameters:
      Double Precision Tolcl1_Umat
!.... Local Variables:
      Integer Istre
      Double Precision Temp1,Temp2
!.... Constants
      Integer Inte1
      Double Precision Real0
      Parameter (Real0=0.0D0,Inte1=1)
      Temp1=Real0
      Temp2=Real0
      Do Istre=Inte1,Nstre
         Temp1=Temp1+(Strea(Istre)-Streb(Istre))**2
         Temp2=Temp1+(Strea(Istre)+Streb(Istre))**2
!Temp2=Temp2+((Strea(Istre)+Streb(Istre))/2+Stora(Istre))**2
      End Do
      If (Abs(Temp2).Le.Real0.Or.Abs(Temp1).Le.Real0) Then
!.... Stress Increment Is Zero
!         Tolcl1=Tolst
         Tolcl1_Umat=Tolst
      Else
         Tolcl1_Umat=Sqrt(Temp1/Temp2)

      End If
      Return
      End

      Function Tolcl2_Umat(Strea,Streb,Nstre,Stora,Tolst)
!....
!     Incremental Stress Norm Function 1
!....
!.... Input Parameters:
      Integer Nstre
      Double Precision Strea(Nstre),Streb(Nstre)
      Double Precision Tolst,Stora(Nstre)

!.... Output Parameters:
      Double Precision Tolcl2_Umat
!.... Local Variables:
      Integer Istre
      Double Precision Temp1,Temp2,Temp3(Nstre),Temp4(Nstre)
      Double Precision Tempr
!.... Constants
      Integer Inte1
      Double Precision Real0,Torr_All,Fbeta
      Parameter (Real0=0.0D0,Inte1=1)
      Temp1=Real0
      Temp2=Real0
      Temp3=Real0
      Temp4=Real0
      Tempr=Real0
      Fbeta=Real0
      Torr_All=0.00000001D0


      Do Istre=Inte1,Nstre
         Temp3(Istre)=(Streb(Istre)-Strea(Istre))*0.5D0
         Temp4(Istre)=Stora(Istre)+(Strea(Istre)+Streb(Istre))*0.5D0
!Temp2=Temp2+((Strea(Istre)+Streb(Istre))/2+Stora(Istre))**2
      End Do


      Do Istre=Inte1,Nstre
         Temp1=Temp1+Temp3(Istre)*Temp3(Istre)
         Temp2=Temp2+Temp4(Istre)*Temp4(Istre)
!Temp2=Temp2+((Strea(Istre)+Streb(Istre))/2+Stora(Istre))**2
      End Do

      Tempr=Sqrt(Temp1/Temp2)

      Tolcl2_Umat=Tempr
      Return
      End
    
    
      Function Tolcl3_Umat(Gausma,Gausmb,Nstre,Lpara,Stora,Tolst,Func)
!....
!     Incremental Stress Norm Function 1
!....
!.... Input Parameters:
      Integer Nstre,Lpara,Func                !Func=1 表示是弹性的判断，Func=1 表示是已知存在塑性应变，而用于计算二三阶仑格库塔的差值
      Double Precision Strea(Nstre),Streb(Nstre),Strec(Nstre)
      Double Precision Parama(Lpara),Paramb(Lpara),Paramc(Lpara)
      Double Precision Gausma(Nstre+Lpara),Gausmb(Nstre+Lpara)
      Double Precision Tolst,Stora(Nstre+Lpara)

!.... Output Parameters:
      Double Precision Tolcl3_Umat
!.... Local Variables:
      Integer Istre
      Double Precision Temp1(Nstre),Temp2(Nstre),Temp3,Temp4
      Double Precision Tempr,Temp5,Temp
!.... Constants
      Integer Inte1
      Double Precision Real0,Torr_All,Fbeta
      Parameter (Real0=0.0D0,Inte1=1)
      Temp1=Real0
      Temp2=Real0
      Temp3=Real0
      Temp4=Real0

      Temp5=Real0

      Tempr=Real0
      Fbeta=Real0
      Torr_All=0.0000000000000001D0

      Goto (1111,1222),Func
!     分割应力和状态变量
      1111  Continue      !如果是1,则用于判断是否进入弹性应力状态
!     分割应力和状态变量
      Do Istre=Inte1,Nstre
         Strea(Istre)=Gausma(Istre)-Stora(Istre)
         Streb(Istre)=Gausmb(Istre)-Gausma(Istre)
         Strec(Istre)=Stora(Istre)
      End Do

      Do Istre=Inte1,Lpara
         Parama(Istre)=Gausma(Nstre+Istre)
         Paramb(Istre)=Gausmb(Nstre+Istre)
         Paramc(Istre)=Stora(Nstre+Istre)
      End Do
!*******************************************************************

      Do Istre=Inte1,Nstre
         Temp1(Istre)=(Streb(Istre)-Strea(Istre))*0.5D0
         Temp2(Istre)=Stora(Istre)+(Strea(Istre)+Streb(Istre))*0.5D0
!Temp2=Temp2+((Strea(Istre)+Streb(Istre))/2+Stora(Istre))**2
      End Do

      Do Istre=Inte1,Nstre
         Temp3=Temp3+Temp1(Istre)*Temp1(Istre)
         Temp4=Temp4+Temp2(Istre)*Temp2(Istre)
      End Do
      Tempr=Sqrt(Temp3/Temp4)
      Temp5=Abs(Parama(2)-Paramb(2))
      If(Tempr.Lt.Torr_All) Tempr=Torr_All     !应力增量的差值
      If(Temp5.Lt.Torr_All) Temp5=Torr_All     !等效塑性应变的差值

      Temp=Max(Tempr,Temp5)
      Tolcl3_Umat=Temp
      Go To 1009

!*******************************************************************
      1222  Continue       !如果是2,则用于计算二阶和三阶龙格库塔的误差
      Temp1=Real0
      Temp2=Real0
      Temp3=Real0
      Temp4=Real0
      Temp5=Real0
      Tempr=Real0
!     分割应力和状态变量
      Do Istre=Inte1,Nstre
         Strea(Istre)=Gausma(Istre)
         Streb(Istre)=Gausmb(Istre)
         Strec(Istre)=Stora(Istre)
      End Do

      Do Istre=Inte1,Lpara
         Parama(Istre)=Gausma(Nstre+Istre)
         Paramb(Istre)=Gausmb(Nstre+Istre)
         Paramc(Istre)=Stora(Nstre+Istre)
      End Do

      Do Istre=Inte1,Nstre
         Temp1(Istre)=(Streb(Istre)-Strea(Istre))*0.5D0
         Temp2(Istre)=Gausmb(Istre)
      End Do

      Do Istre=Inte1,Nstre
         Temp3=Temp3+Temp1(Istre)*Temp1(Istre)
         Temp4=Temp4+Temp2(Istre)*Temp2(Istre)
      End Do

      Tempr=Sqrt(Temp3/Temp4)
      Temp5=Abs(Parama(2)-Paramb(2))
      If(Tempr.Lt.Torr_All) Tempr=Torr_All     !应力增量的差值
      If(Temp5.Lt.Torr_All) Temp5=Torr_All     !等效塑性应变的差值
      Tolcl3_Umat=Max(Tempr,Temp5)

      Go To 1009
!*******************************************************************
      1009  Continue
      Return
      End
    
    
      Logical Function D0Eq(A,B)
      Implicit None
      Double Precision A,B,Err,Dmach_Umat
      External Dmach_Umat
      Logical D0Eq0_Umat
      Save Err
      Data Err/-1.0D0/
      If (Err.Le.0.0D0) Err=100.0D0*Dmach_Umat(1)
      D0Eq=Abs(A-B).Le.Err*Abs(A+B).Or.(D0Eq0_Umat(A).And.D0Eq0_Umat(B))
      Return
      End

      Logical Function D0Eq0_Umat(A)
      Implicit None
      !      Include 'Kindfile'
      Double Precision A,Err,Dmach_Umat
      External Dmach_Umat
      Save Err
      Data Err/-1.0D0/
      If (Err.Le.0.0D0) Err=100.0D0*Dmach_Umat(2)
      D0Eq0_Umat=Abs(A).Le.Err
      Return
      End

      Function Dmach_Umat(Job)
!      Include 'Kindfile'
      Double Precision Dmach_Umat
      Integer Job
!
!     Smach Computes Machine Parameters Of Floating Point
!
      Double Precision Eps,Tiny,Huge,S
!
      Eps = 1.0D0
      10 Eps = Eps/2.0D0
      S = 1.0D0 + Eps
      Call Dummy_Umat(Eps)
      If (S .Gt. 1.0D0) Go To 10
      Eps = 2.0D0*Eps
!
      S = 1.0D0
      20 Tiny = S
      S = S/16.0D0
      Call Dummy_Umat(S)
      If (Abs(S*1.0) .Gt. 0.0D0) Go To 20
      Tiny = (Tiny/Eps)*100.0
      Huge = 1.0D0/Tiny
!
      If (Job .Eq. 1) Dmach_Umat = Eps
      If (Job .Eq. 2) Dmach_Umat = Tiny
      If (Job .Eq. 3) Dmach_Umat = Huge
      Return
      End


      Subroutine Dummy_Umat(R)
!.... Dummy_Umat Subroutine To Avoid Optimization Of Key Variable
!      Include 'Kindfile'
      Double Precision R,S
      Common /Dummy_Umat1/ S
      S=R
      Return
      End


      Subroutine Dcopy_Umat(N,Dx,Incx,Dy,Incy)
!
!  -- Reference Blas Level1 Routine --
!  -- Reference Blas Is A Software Package Provided By Univ. Of Tennessee,    --
!  -- Univ. Of California Berkeley, Univ. Of Colorado Denver And Nag Ltd..--
!
!     .. Scalar Arguments ..
      Integer Incx,Incy,N
!     ..
!     .. Array Arguments ..
      Double Precision Dx(*),Dy(*)
!
!  =====================================================================
!     .. Local Scalars ..
      Integer I,Ix,Iy,M,Mp1
!     ..
!     .. Intrinsic Functions ..
      Intrinsic Mod
!     ..
      If (N.Le.0) Return
      If (Incx.Eq.1 .And. Incy.Eq.1) Then
!
!        Code For Both Increments Equal To 1
!        Clean-Up Loop
!
         M = Mod(N,7)
         If (M.Ne.0) Then
            Do I = 1,M
               Dy(I) = Dx(I)
            End Do
            If (N.Lt.7) Return
         End If
         Mp1 = M + 1
         Do I = Mp1,N,7
            Dy(I) = Dx(I)
            Dy(I+1) = Dx(I+1)
            Dy(I+2) = Dx(I+2)
            Dy(I+3) = Dx(I+3)
            Dy(I+4) = Dx(I+4)
            Dy(I+5) = Dx(I+5)
            Dy(I+6) = Dx(I+6)
         End Do
      Else
!
!        Code For Unequal Increments Or Equal Increments
!          Not Equal To 1
!
         Ix = 1
         Iy = 1
         If (Incx.Lt.0) Ix = (-N+1)*Incx + 1
         If (Incy.Lt.0) Iy = (-N+1)*Incy + 1
         Do I = 1,N
            Dy(Iy) = Dx(Ix)
            Ix = Ix + Incx
            Iy = Iy + Incy
         End Do
      End If
      Return
!
!     End Of Dcopy_Umat
!
      End

!
!!********************************************************************************
!!            *************** Get Stress Finished*************************
!!********************************************************************************
!
!
!********************************************************************************
!            *************** Get Dematrix *************************
!********************************************************************************

      Subroutine Getdmatrx(Propd,Lprpd,Param,Lpara,Stran,Dstan,Stress,Dstre,Dmatx,Nstre,Ielem,Igaus,Ktest,Ntest)

      Implicit None
!     Associative Classical Model With Cohesion And Variable
!     Elastic Modulus, Elasto Plasticity
!     Taken From Owen And Hinton - Finite Element In Plasticity

!     Input Variables
      Integer Lprpd,Iswdp,Icdat,Icout,Lpara,Ielem,Igaus,Ntest
      Integer Nstre
      Integer Ktest(Ntest)

      Double Precision Voidr,Stran(Nstre)

!     Input/Output Variables

      Double Precision Propd(Lprpd)
      Double Precision Param(Lpara)

      Double Precision Indstan(Nstre),Indstre(Nstre),Stress(Nstre)
      Double Precision Dstan(Nstre),Dstre(Nstre)
      Double Precision Estre(Nstre),Dmatx(Nstre,Nstre)
      Integer::I,J

      Icdat=1
      Icout=1

      Iswdp=4
      Call Mcoul3_Umat(Propd,Lprpd,Iswdp,Icdat,Icout,Dstre,Dstan,Dmatx,Stress,   &
                       Param,Lpara,Ielem,Igaus,Voidr,Ktest,Ntest,Nstre)

      Return
      End Subroutine


!********************************************************************************
!            *************** Get Fssi Mc *************************
!********************************************************************************
      Subroutine Mcoul3_Umat(Propd,Lprpd,Iswdp,Icdat,Icout,Indstre,Indstan,Inoutdmatx,   &
                             Inestre,Param,Lpara,Ielem,Igaus,Voidr,Ktest,Ntest,Nstre)
      Implicit None

!     Associative Classical Model With Cohesion And Variable
!     Elastic Modulus, Elasto Plasticity
!     Taken From Owen And Hinton - Finite Element In Plasticity

!     Input Variables
      Integer Lprpd,Iswdp,Icdat,Icout,Lpara,Ielem,Igaus,Ntest
      Integer Nstre
      Integer Ktest(Ntest)

      Double Precision Voidr,Dflag

!     Input/Output Variables

      Double Precision Propd(Lprpd)
      Double Precision Param(Lpara)
      Double Precision Dmatx_Inv(6,6)      ! 弹性D矩阵的逆矩阵
      Double Precision PlasticStrain(6)    ! 塑性应变

!**************************************************************
!  类    型：更改代码
!  功    能：将传进来的 Dstre，Dstan，Dstan 改为新的数组，目的是
!            需要根据维度来对上述三个数组进行扩充，为避免对后续
!            代码进行大的改正，因此重新定义传进来的量。分别命名为：
!            Indstre，Indstan，Indstan  意义为：传进来的数据
!  源 代 码:
!            Double Precision Dstan(Nstre),Dstre(Nstre),Estre(Nstre)
!  更改代码：
      Double Precision Indstan(Nstre),Indstre(Nstre),Inestre(Nstre),Inoutdmatx(Nstre,Nstre)
      Double Precision Dstan(6),Dstre(6),Estre(6),Dmatx(6,6)
      Integer::I,J
!  备    注：   王天民20211123
!**************************************************************

!     Propd Array
!     1: Young'S Modulus
!     2: Poisson'S Ratio
!     3: Cohesion
!     4: Work Hardening Modulus As Defined In The Book
!     5: Friction Angle
!     6: Criterion Number: 1. Tresca
!                          2. Von Mises
!                          3. Mohr Coulomb
!                          4. Drucker Prager
!                          5. 采用陈卫忠老师书上给出的Mc 模型,对传统Mc 的棱角进行光滑处理，但是积分算法仍然是显示积分算法
!     7: 选择积分算法

!     Param Array:
!     1. Effst
!     2. Epstn
!     3. Escur

!     Local Variables
      Integer Istre,Ncrit
      Double Precision Fstre(Nstre),Fstre1(6)
      Integer Kelem,Lelem,Kgaus,Lgaus,Kparm,Lparm,Nbegn,Nlast,Nfunc


      Kelem=-1
      Lelem=-1
      Kgaus=-1
      Lgaus=-1
      Kparm=-1
      Lparm=-1
      Nbegn=-1
      Nlast=-1
      Nfunc=-1
      Dflag=0.0D0

!**************************************************************
!  新增代码：
!  功    能：将传进来的总应变，增量应变，总应力根据维数进行扩充
!            当Nstre=4 时，即平面应变问题时，本构程序中依然按照
!            3维进行处理，只需将后面的两个应力补位为0 即可；
!            注意：Fssicas 中的应力顺序为：Sx Sy Sz Txy Tyz Txz
      Dstan=0.0D0
      Dstre=0.0D0
      Estre=0.0D0
      Dmatx=0.0D0
      If (Int(Nstre).Eq.6)Then
         Dstan=Indstan
         Dstre=Indstre
         Estre=Inestre
      Else If (Int(Nstre).Eq.4)Then
         Do I=1,4
            Dstan(I)=Indstan(I)
            Dstre(I)=Indstre(I)
            Estre(I)=Inestre(I)
         End Do
      Else
         Print *,'Plane Stress Problem Has Not Been Considered'
      End If

!  备    注：   王天民20211123
!**************************************************************

!     Go To The Choice Of Action
      Goto (1100,1200,1300,1400,1500),Iswdp

1100  Continue

1200  Continue
!     Calculate The Incremental Stress
!     Save The Original Stress State

      Fstre1=Estre

      Ncrit=Nint(Propd(6))

      Call Plast2_Umat(Dstan,Estre,Icout,Propd(1),Propd(2),Propd(3),Propd(4),Propd(5),Param(1),Param(2),Param(3),Ncrit)

      Do Istre=1,6
         Dstre(Istre)=Estre(Istre)-Fstre1(Istre)
      End Do

!     求塑性应变：PlasticStrain=增量应变-增量应力/弹性D矩阵 ，用于在状态变量中输出塑性总应变 [何坤鹏 2024年3月]
      Call Modps1_Umat(Propd(1),Propd(2),Dmatx)
      Call Matrix_Inverted_Mkl(Dmatx,Dmatx_Inv,6,1)
      PlasticStrain=Dstan-Matmul(Dstre,Dmatx_Inv)
      
      Param(3+1:3+Nstre)=Param(3+1:3+Nstre)+PlasticStrain
      
!**************************************************************
!  功  能： 根据维度，输出对应维度的增量应力和总应力
      If (Int(Nstre).Eq.6)Then
         Indstre=Dstre
         Inestre=Estre
      Else If (Int(Nstre).Eq.4)Then
         Do I=1,4
            Indstre(I)=Dstre(I)
            Inestre(I)=Estre(I)
         End Do
      Else
         Print *,'There Is A Problem At The Progress Of Form Stress And Dstress'
      End If
      Goto 1900
!  备  注：   王天民20211123
!**************************************************************


1300 Continue
!  Form Symmetric D-Matrix
!**************************************************************
!  功  能：   将总应力 Szz 和 Txy 调换位置
      Ncrit=Nint(Propd(6))
      Call Plast3_Umat(Dmatx,Estre,Icout,Propd(1),Propd(2),Propd(3),Propd(4),Propd(5),Param(1),Param(2),Param(3),Ncrit,Dflag)
      Goto 1450

!  备  注：   王天民20211123
!**************************************************************


1400  Continue
!  Form Not Necessary Symmetric D-Matrix
!**************************************************************
!  功  能：   将总应力 Szz 和 Txy 调换位置
      Ncrit=Nint(Propd(6))
      Call Plast3_Umat(Dmatx,Estre,Icout,Propd(1),Propd(2),Propd(3),Propd(4),Propd(5),Param(1),Param(2),Param(3),Ncrit,Dflag)

!  备  注：   王天民20211123
!**************************************************************

1450  Continue
!**************************************************************
!  功  能： 根据维度，输出对应维度的刚度矩阵

      If (Int(Nstre).Eq.6)Then
         Indstre=Dstre
         Inestre=Estre
      Else If (Int(Nstre).Eq.4)Then
         Do I=1,4
            Indstre(I)=Dstre(I)
            Inestre(I)=Estre(I)
         End Do
      Else
         Print *,'There Is A Problem At The Progress Of Form Stress And Dstress'
      End If

      If (Int(Nstre).Eq.6)Then
         Inoutdmatx=Dmatx
      Else If (Int(Nstre).Eq.4)Then
         Do I=1,4
            Do J=1,4
               Inoutdmatx(I,J)=Dmatx(I,J)
            End Do
         End Do
      Else
         Print *,'There Is A Problem At The Progress Of Form Dmatx'
      End If
      Goto 1900
!   备  注：   王天民20211123
!**************************************************************

1500 Continue   ! Initialization
!     Create A Zero Incremental Strain So That Plast2 Can Be Used
!     To Calculate Escur
!     Do Istre=1,Nstre
!        Dstan(Istre)=0.0D0
!     End Do
!     Param(1)=0.0D0
!     Param(2)=0.0D0

      Goto 1900

1900  Continue   ! Common Exit

      Return
      End


      Subroutine Plast2_Umat(Stran,Strsg,Icout,Young,Poiss,Cohes,Hards,Frict,Effst,Epstn,Escur,Ncrit)
!     Update Stress State Due To Strain Increment
      Implicit None

!     Incremental Strain
      Double Precision Stran
!     Current Stress State
      Double Precision Strsg
!     Output Channel
      Integer Icout
!     Young'S Modulus
      Double Precision Young
!     Poisson'S Ratio
      Double Precision Poiss
!     Uniaxial Yield Stress
      Double Precision Uniax
!     One Dimensional Yield Stress
      Double Precision Cohes
!     Hardening Modulus
      Double Precision Hards
!     Friction Angle In Degree
      Double Precision Frict
!     Effective Stress State
      Double Precision Effst
!     Elastoplastic Strain
      Double Precision Epstn
!     Current Stress State
      Double Precision Escur
!     Criteria Number
      Integer Ncrit
      Double Precision Real0,Real1,Real3,Real6,Real8
      Integer Nstre,Nstr1
      Integer Kunld
!     Elasto-Plastic Matrix
      Double Precision Dmatx
      Integer Istre
!     Incremental Stress
      Double Precision Stres
      Integer Jstre
!     Previous Yield Stress
      Double Precision Preys
      Integer Istr1
!     Incremental Elastic Stress
      Double Precision Desig
!     Deviatoric Stress
      Double Precision Devia
!     A-Vector
      Double Precision Avect
!     D-Vector
      Double Precision Dvect
!     Current Total Stress
      Double Precision Sgtot
!     Current Stress State
      Double Precision Sigma
      Double Precision Sint3,Steff,Theta,Varj2
      Double Precision Yield,Espre,Rfact,Astep,Reduc
      Integer Mstep,Istep
      Double Precision Abeta,Agash,Dlamd,Bgash,Curys
      Double Precision Bring
      Double Precision Root3
      Parameter(Real0=0.0D0, Real1=1.0D0, Real3=3.0D0,Real6=6.0D0,Real8=8.0D0)

!**************************************************************
!  新增代码：
!  功    能：将传进来的总应变，增量应变，总应力根据维数进行扩充
!  源 代 码：Parameter(Nstre=4,Nstr1=4)
      Parameter(Nstre=6,Nstr1=6)
      Dimension Stran(Nstre),Strsg(Nstre),Dmatx(Nstre,Nstre),Desig(Nstre),Sigma(Nstre),   &
                Devia(Nstre),Avect(Nstre),Dvect(Nstre),Sgtot(Nstre),Stres(Nstre)

!  备    注：   王天民20211123
!**************************************************************

      Abeta=0.0D0
      Root3=Sqrt(3.0D0)
      Kunld=0
      Uniax=Cohes

      If(Ncrit.Eq.3)Uniax=Cohes*Cos(Frict*Atan(1.0D0)/45.0D0)
      If(Ncrit.Eq.4)Uniax=Real6*Cohes*Cos(Frict*Atan(1.0D0)/45.0D0)/(Root3*(Real3-Sin(Frict*Atan(1.0D0)/45.0D0)))
      If(Ncrit.Eq.5)Uniax=Cohes*Cos(Frict*Atan(1.0D0)/45.0D0)

!     Set Up Elastic D Matrix And Calculate Stress Assuming Linearity
      Call Modps1_Umat(Young,Poiss,Dmatx)
      Do Istre=1,Nstre
         Stres(Istre)=Real0
         Do Jstre=1,Nstre
            Stres(Istre)=Stres(Istre)+Dmatx(Istre,Jstre)*Stran(Jstre)
         End Do
      End Do

      Preys=Uniax+Epstn*Hards
      Do Istr1=1,Nstr1
         Desig(Istr1)=Stres(Istr1)
         Sigma(Istr1)=Strsg(Istr1)+Stres(Istr1)
      End Do
      Call Invariants_Mohrcoulomb_Umat(Devia,Frict,Ncrit,Sint3,Steff,Sigma,Theta,Varj2,Yield,Nstre)
      Espre = Effst - Preys             !检查该高斯点上的应力在上一步迭代中是否进入屈服
      If(Espre.Ge.Real0)Go To 30        !屈服值和当前函数值的差异
      Escur = Yield - Preys
      If(Escur.Le.Real0)Go To 40
      Rfact = Escur/(Yield-Effst)
      Go To 50
 30   Escur = Yield - Effst
      If(Escur.Gt.Real0)Go To 35
      Go To 40
 35   Rfact = Real1

 50   Astep = Escur * Real8/Uniax+Real1
      If (Astep.Gt.20) Then
         Mstep=100
      Else
         Mstep=Int(Astep)
      End If

      Astep = Mstep
      Reduc = Real1 - Rfact
      Do Istr1 = 1,Nstr1
         Sgtot(Istr1)=Strsg(Istr1)+Reduc*Stres(Istr1)
         Stres(Istr1)=Rfact*Stres(Istr1)/Astep
      End Do

!     将应力增量分成很多份，逐份拉回
      Do Istep=1,Mstep
         Call Invariants_Mohrcoulomb_Umat(Devia,Frict,Ncrit,Sint3,Steff,Sgtot,Theta,Varj2,Yield,Nstre) ! Obtain-> Varj2,Yield
         Call Yieldf_Umat(Avect,Devia,Frict,Ncrit,Nstr1,Steff,Theta,Varj2)                             ! Evaluation Of Flow Vectors->Avect 
         Call Flowpl_Umat(Avect,Abeta,Dvect,Nstr1,Young,Poiss,Hards)                                   ! Obtain-> Abeta,Dvect
         Agash=Real0
         Do Istr1=1,Nstr1
            Agash = Agash + Avect(Istr1)*Stres(Istr1)
         End Do
         Dlamd = Agash*Abeta
         If(Dlamd.Lt.Real0)Dlamd=Real0        ! 说明是中性变载，那么塑性模量等于0
         Bgash=Real0
         Do Istre=1,Nstre
            Bgash = Bgash + Avect(Istre)*Sgtot(Istre)
            Sgtot(Istre)=Sgtot(Istre) + Stres(Istre)- Dlamd*Dvect(Istre)
         End Do
         
         !WTM20240614 解决C值为0时，不能计算的问题
         If(Abs(Yield).Lt.1.0d-10) Yield = 1.0d-10 
         
         Epstn = Epstn + Dlamd*Bgash/Yield
      End Do

      Call Invariants_Mohrcoulomb_Umat(Devia,Frict,Ncrit,Sint3,Steff,Sgtot,Theta,Varj2,Yield,Nstre)
      Curys = Uniax + Epstn*Hards
      Bring =1.0D0
      
      !WTM20240614 解决C值为0时，不能计算的问题
      If(Abs(Yield).Lt.1.0d-10) Yield = 1.0d-10
      
      If(Yield.Gt.Curys) Bring = Curys/Yield
      Do Istr1=1,Nstr1
         Sgtot(Istr1)=Bring*Sgtot(Istr1)
      End Do
      Effst = Bring*Yield
      Strsg = Sgtot

!     End Of Inclusion
      Go To 190

 40   Do Istr1=1,Nstr1
         Strsg(Istr1)=Strsg(Istr1)+Desig(Istr1)
      End Do
      Effst = Yield

190   Continue

      Return
      End


      Subroutine Plast3_Umat(Dmatx,Strsg,Icout,Young,Poiss,Cohes,Hards,Frict,Effst,Epstn,Escur,Ncrit,Dflag)
!     To Form Elastoplastic D Matrix
      Implicit None
      Double Precision Dmatx,Dflag
      Double Precision Strsg
      Integer Icout
      Double Precision Young
      Double Precision Poiss
      Double Precision Cohes
      Double Precision Hards
      Double Precision Frict
      Double Precision Epstn
      Double Precision Escur
      Integer Ncrit
      Integer Nstre
      Integer Nstr1
      Double Precision Devia
      Double Precision Avect
      Double Precision Dvect
      Double Precision Stres
!**************************************************************
!  源 代 码：Parameter(Nstre=4,Nstr1=4)
      Parameter(Nstre=6,Nstr1=6)
!  备    注：   王天民20211123
!**************************************************************
      Double Precision Real0,Real1
      Integer Istr1
      Double Precision Sint3,Steff,Theta,Varj2,Yield,Abeta
      Integer Istre,Jstre
      Double Precision Effst
!**************************************************************
!  源 代 码：
      Dimension Dmatx(Nstre,Nstre),Strsg(Nstre),Stres(Nstre),Dvect(Nstre),Avect(Nstre),Devia(Nstre)
!      Double Precision Dmatx(Nstre,Nstre),Strsg(Nstre),Stres(Nstre),Dvect(Nstre),Avect(Nstre),Devia(Nstre)
!  备    注：   王天民20211123
!**************************************************************

      Parameter(Real0=0.0D0,Real1=1.0D0)
      Double Precision Vsmal
      Parameter (Vsmal=1.0D-30)

      Call Modps1_Umat(Young,Poiss,Dmatx)
      If(Escur.Le.Real0.Or.Abs(Dflag-Real1).Le.Vsmal)Go To 10
      Do Istr1=1,Nstr1
         Stres(Istr1) = Strsg(Istr1)
      End Do

      Call Invariants_Mohrcoulomb_Umat(Devia,Frict,Ncrit,Sint3,Steff,Stres,Theta,Varj2,Yield,Nstre)
      Call Yieldf_Umat(Avect,Devia,Frict,Ncrit,Nstr1,Steff,Theta,Varj2)
      Call Flowpl_Umat(Avect,Abeta,Dvect,Nstr1,Young,Poiss,Hards)

      Do Istre = 1,Nstre
         Do Jstre=1,Nstre
            Dmatx(Istre,Jstre)=Dmatx(Istre,Jstre) - Abeta*Dvect(Istre)*Dvect(Jstre)
         End Do
      End Do

      10   Continue

      Return
      End


      Subroutine Flowpl_Umat(Avect,Abeta,Dvect,Nstr1,Young,Poiss,Hards)
!     Multiplies Avect With Dmatx, And Obtains Abeta
      Implicit None
      Double Precision Avect
      Double Precision Abeta
      Double Precision Dvect
      Integer Nstr1
      Integer Nstre
      Double Precision Real0
      Double Precision Real1,Real2
      Double Precision Young
      Double Precision Poiss
      Double Precision Hards
      Double Precision Fmul1
      Double Precision Fmul2
      Double Precision Denom
      Integer Istr1
!**************************************************************
!  源 代 码：Parameter(Nstre=4)
      Parameter(Nstre=6)                        !  备    注：   王天民20211123
!**************************************************************
      Parameter(Real0=0.0D0)
      Parameter (Real1=1.0D0,Real2=2.0D0)
!**************************************************************
!  源 代 码：
      Dimension Avect(Nstre), Dvect(Nstre)
!      Double Precision Avect(Nstre), Dvect(Nstre)!  备    注：   王天民20211123
!**************************************************************
      Fmul1=Young/(Real1+Poiss)
      Fmul2=Young*Poiss*(Avect(1)+Avect(2)+Avect(3))/((Real1+Poiss)*(Real1-Real2*Poiss))
      Dvect(1)=Fmul1*Avect(1)+Fmul2
      Dvect(2)=Fmul1*Avect(2)+Fmul2
      Dvect(3)=Fmul1*Avect(3)+Fmul2
      Dvect(4)=0.5D0*Avect(4)*Young/(Real1+Poiss)
      Dvect(5)=0.5D0*Avect(5)*Young/(Real1+Poiss)
      Dvect(6)=0.5D0*Avect(6)*Young/(Real1+Poiss)

      Denom = Hards
      Do Istr1=1,Nstre
         Denom = Denom + Avect(Istr1)*Dvect(Istr1)
      End Do
      Abeta = Real1/Denom

      Return
      End


      Subroutine Modps1_Umat(Young,Poiss,Dmatx)

!     Sets Up Elastic D Matrix
      Implicit None
      Integer Istr1
      Integer Jstr1
      Integer Nstre
      Double Precision Real0
      Double Precision Real1
      Double Precision Real2
      Double Precision Young
      Double Precision Poiss
      Double Precision Dmatx
      Double Precision Const
      Double Precision Conss
      Parameter(Real0=0.0D0,Real1=1.0D0,Real2=2.0D0)
!**************************************************************
!  源 代 码：Parameter(Nstre=4)
      Parameter(Nstre=6)
      Dimension Dmatx(Nstre,Nstre)
!  备    注：   王天民20211123
!      Double Precision Dmatx(Nstre,Nstre)
!**************************************************************


!      Save    ! 影响并行【1.24】

!     Initialise The Array
      Do Istr1=1,Nstre
         Do Jstr1=1,Nstre
            Dmatx(Istr1,Jstr1)=Real0
         End Do
      End Do
!
      Const = Young*(Real1-Poiss)/((Real1+Poiss)*(Real1-Real2*Poiss))
      Conss= Poiss/(Real1-Poiss)
      Dmatx(1,1)=Const
      Dmatx(2,2)=Const
      Dmatx(3,3)=Const
      Dmatx(1,2)=Const*Conss
      Dmatx(2,1)=Const*Conss
      Dmatx(1,3)=Const*Conss
      Dmatx(3,1)=Const*Conss
      Dmatx(2,3)=Const*Conss
      Dmatx(3,2)=Const*Conss
      Dmatx(4,4)=Const*(Real1-Real2*Poiss)/(Real2*(Real1-Poiss))
      Dmatx(5,5)=Const*(Real1-Real2*Poiss)/(Real2*(Real1-Poiss))
      Dmatx(6,6)=Const*(Real1-Real2*Poiss)/(Real2*(Real1-Poiss))

      Return
      End


      Subroutine Invariants_Mohrcoulomb_Umat(Devia,Frict,Ncrit,Sint3,Steff,Stemp,Theta,Varj2,Yield,Nstre)
!     Evaluates Stress Invariant And Value Of The Yield Function

      Implicit None

      Integer Ncrit,I,Nstre
      Double Precision Devia(Nstre)
      Double Precision Frict

      Double Precision Sint3
      Double Precision Steff
      Double Precision Stemp(Nstre)
      Double Precision Theta
      Double Precision Varj2
      Double Precision Yield
      Double Precision Root3
      Double Precision Smean
      Double Precision Varj3
      Double Precision Phira
      Double Precision Snphi
      Double Precision Real0,Real1,Real2,Real6,Real3


!**************************************************************
!  新增代码：
!  功    能：将全局变量重新定义，避免后面数据交换
!  源 代 码：
!      Double Precision Devia(Nstre),Stemp(Nstre)
!  备    注：   王天民20211123
!**************************************************************

      Double Precision Vsmal
      Parameter (Vsmal=1.0D-30)
      Double Precision Sterm
      Double Precision Acon1,Acon2    ! 这里的 Acon1,Acon2 分别表示 ，平均应力Sm前的系数、Sqrt(J2)前的系数
      Double Precision::Thetat,A1,A2,B1,B2,Kf,Sgn
      Double Precision::Const1,Const2,S3Ta,Sbar
      Double Precision One,Two,Three,Pi
      Parameter(One=1.0D0,Two=2.0D0,Three=3.0D0)


      Real0=0.0D0
      Real1=1.0D0
      Real2=2.0D0
      Real6=6.0D0
      Real3=3.0D0
      Pi=3.1415927

      Root3=Sqrt(3.0D0)
!**************************************************************
!  功    能：将传进来的总应变，增量应变，总应力根据维数进行扩充
      Smean=(Stemp(1)+Stemp(2)+Stemp(3))/Real3
      Devia(1)=Stemp(1)-Smean     !Sx
      Devia(2)=Stemp(2)-Smean     !Sy
      Devia(3)=Stemp(3)-Smean     !Sz
      Devia(4)=Stemp(4)           !Txy
      Devia(5)=Stemp(5)           !Tyz
      Devia(6)=Stemp(6)           !Tzx
      Varj2=0.5*(Devia(1)*Devia(1)+Devia(2)*Devia(2)+Devia(3)*Devia(3))+Devia(4)*Devia(4)+Devia(5)*Devia(5)+Devia(6)*Devia(6)
      Varj3=Devia(1)*Devia(2)*Devia(3)+2*Devia(4)*Devia(5)*Devia(6)-Devia(1)*Devia(5)*Devia(5)-Devia(2)*Devia(6)*Devia(6)-Devia(3)*Devia(4)*Devia(4)
!  备    注：   王天民20211123
!**************************************************************
!
!If (Ieee_Is_Nan(Varj2)) Then
!   Print *,'Stress:',Stemp
!   Print *,'Varj2 In Invar:',Varj2
!   Stop 'Stopped In Invar Due To Varj2 Not A Number'
!End If

      Steff=Sqrt(Varj2)
      If(Abs(Steff).Le.Vsmal)Go To 10
      If (Abs(Real2*Varj2*Steff).Gt.0.0D0) Then
         Sint3=-Real3*Root3*Varj3/(Real2*Varj2*Steff)
      Else
         Sint3=Real1
      End If
      If(Sint3.Gt.Real1)Sint3=Real1
      Go To 20
      10        Sint3=Real0

      20        Continue

      If(Sint3.Lt.-Real1)Sint3=-Real1
      If(Sint3.Gt.Real1)Sint3=Real1
      Theta=Asin(Sint3)/Real3
      Go To (1,2,3,4,5) Ncrit
      Stop 'Stopped In Invar Due To Ncrit'

!     Tresca
      1     Acon1=0.0D0
      Acon2=Real2*Cos(Theta)
      Yield=Acon2*Steff
      Return



!     Von Mises
      2     Acon1=0.0D0
      Acon2=Root3
      Sterm=Steff*Acon2
      Yield=Acon2*Steff
      Return

!     Mohr Coulomb
      3     Phira=Frict*Atan(1.0D0)/45.0D0
      Snphi=Sin(Phira)
      Acon1=Snphi
      Acon2=Cos(Theta)-Sin(Theta)*Snphi/Root3
      Sterm=Steff*Acon2
      Yield=Acon1*Smean+Acon2*Steff
      Return

!     Drucker Prager
!     外接圆锥（外角点外接圆）！ 详见郑颖人《岩土塑性力学第二版》P74页 表4-1 Dp1
      4     Phira=Frict*Atan(1.0D0)/45.0D0
      Snphi=Sin(Phira)
      Acon1=Real6*Snphi/(Root3*(Real3-Snphi))
      Acon2=1.0D0
      Sterm=Steff*Acon2
      Yield=Acon1*Smean+Acon2*Steff        !注意此处是Smean  不是I1
      Return


!     Mohr Coulomb Chen

      5     Phira=Frict*Atan(1.0D0)/45.0D0
      Snphi=Sin(Phira)
      Acon1=Snphi

!!    Chen_Mc 的屈服函数
      S3Ta=Sin(Three*Theta)
      Theta=Asin(S3Ta)/Three
      Thetat=(25.0/180.0)*Pi
      A1=1.432052062044D0
      A2=0.406941858374D0
      B1=0.544290524902D0
      B2=0.673903324498D0
      If (Abs(Theta).Le.Thetat) Then
         Kf=Cos(Theta)-Sin(Phira)*Sin(Theta)/Sqrt(Three)
      Else
         Sgn=Sign(One,Theta)
         Const1=A1+A2*Sgn*Sin(Phira)
         Const2=B1*Sgn+B2*Sin(Phira)
         Kf=Const1-Const2*S3Ta
      End If

      Sbar=Steff
      Yield=Acon1*Smean+Sbar*Kf
      Return

      End
!

      Subroutine Yieldf_Umat(Avect,Devia,Frict,Ncrit,Nstr1,Steff,Theta,Varj2)
!     Evaluation Of Flow Vectors->Avect

      Implicit None
      Double Precision Avect
      Double Precision Devia
      Double Precision Frict
      Integer Ncrit
      Integer Nstr1
      Double Precision Steff
      Integer Nstre
      Double Precision Theta
      Double Precision Varj2
      Double Precision Real1,Real0,Real2,Real3
      Double Precision Three
      Double Precision Veca1
      Double Precision Veca2
      Double Precision Veca3
      Double Precision Tanth
      Double Precision Tant3
      Double Precision Sinth
      Double Precision Costh
      Double Precision Cost3
      Double Precision Root3
      Integer Istr1
      Double Precision Cons1
      Double Precision Abthe
      Double Precision Cons2
      Double Precision Cons3
      Double Precision Plumi
      Double Precision Snphi
      Double Precision Sinfri
      Logical Isotr

      Double Precision One,Two,Pi
!**************************************************************
!  类    型：更改代码
!  功    能：扩展到三维
!  源 代 码：Parameter(Nstre=4)
      Parameter(Nstre=6)

!  备    注：   王天民20211123
!**************************************************************

      Parameter(Real1=1.0D0,Real0=0.0D0,Real2=2.0D0,Real3=3.0D0)
      Parameter(One=1.0D0,Two=2.0D0,Three=3.0D0,Pi=3.1415927)
!**************************************************************
!  类    型：更改代码
!  功    能：将全局变量进行更改
!  源 代 码：
      Dimension Avect(Nstre),Devia(Nstre),Veca1(Nstre),Veca2(Nstre),Veca3(Nstre)
!     Double Precision Avect(Nstre),Devia(Nstre),Veca1(Nstre),Veca2(Nstre),Veca3(Nstre)
!  备    注：   王天民20211123
!**************************************************************
      Double Precision::Thetat,A1,A2,B1,B2,Kf,Sgn,Sbar,Zero
      Double Precision::Const1,Const2,S3Ta,Dkf,Alfa,Phira

      Double Precision Vsmal
      Parameter (Vsmal=1.0D-30)
!
      Isotr=Abs(Steff).Le.Vsmal
      Tanth=Tan(Theta)
      Tant3=Tan(Real3*Theta)
      Sinth=Sin(Theta)
      Costh=Cos(Theta)
      Cost3=Cos(Real3*Theta)
      Root3=Sqrt(3.0D0)
      Zero=0.0D0


!**************************************************************
!  功    能：扩展到三维
!     Calculate Vector A1
      Veca1(1)=Real1
      Veca1(2)=Real1
      Veca1(3)=Real1
      Veca1(4)=Real0
      Veca1(5)=Real0
      Veca1(6)=Real0

!     Calculate Vector A2
      If (Isotr)Then
         Veca2=0.0D0
      Else
         Do Istr1=1,Nstr1
            Veca2(Istr1)=Devia(Istr1)/(Real2*Steff)
         End Do
         Veca2(4)=Devia(4)/Steff
         Veca2(5)=Devia(5)/Steff
         Veca2(6)=Devia(6)/Steff
      End If
!     Calculate Vector A3

      Veca3(1)=Devia(2)*Devia(3)-Devia(5)*Devia(5)+Varj2/Real3
      Veca3(2)=Devia(3)*Devia(1)-Devia(6)*Devia(6)+Varj2/Real3
      Veca3(3)=Devia(1)*Devia(2)-Devia(4)*Devia(4)+Varj2/Real3
      Veca3(4)=Real2*(Devia(5)*Devia(6)-Devia(3)*Devia(4))
      Veca3(5)=Real2*(Devia(6)*Devia(4)-Devia(1)*Devia(5))         ! 根据 Pz3 程序中，后面三项是 Txy Tyz Tzx
      Veca3(6)=Real2*(Devia(4)*Devia(5)-Devia(2)*Devia(6))
!  备    注：   王天民20211123
!**************************************************************

      Goto (1,2,3,4,5) Ncrit

!** Tresca
!
      1  Cons1=Real0
      Abthe=Abs(Theta*57.29577951308D0)
      If(Abthe.Lt.29.999999D0)Go To 20
      Cons2=Root3
      Cons3=Real0
      Go To 40
      20  Cons2=Real2*(Costh+Sinth*Tant3)
      Cons3=Root3*Sinth/(Varj2*Cost3)
      Go To 40
!
!** Von Mises
!
      2  Cons1=Real0
      Cons2=Root3
      Cons3=Real0
      Go To 40
!
!** Mohr Coulomb
!
      3   Cons1=Sin(Frict*Atan(1.0D0)/45.0D0)/Real3
      Abthe=Abs(Theta*57.29577951308D0)
      If(Abthe.Lt.29.999999D0)Go To 30
      Cons2=Real0
      Cons3=Real0
      Plumi=Real1
      If(Theta.Gt.Real0)Plumi=-Real1
      Cons2=0.5D0*(Root3+Plumi*Cons1*Root3)
      Go To 40
      30  Continue
      If (Isotr) Then
         Cons2=Real0
         Cons3=Real0
      Else
         Cons2=Costh*((Real1+Tanth*Tant3)+Cons1*(Tant3-Tanth)*Root3)
         Cons3=(Root3*Sinth+Real3*Cons1*Costh)/(Real2*Varj2*Cost3)
      End If
      Go To 40
!
!**  Drucker Prager
!
      4   Snphi=Sin(Frict*Atan(1.0D0)/45.0D0)
      Cons1=Real2*Snphi/(Root3*(Real3-Snphi))
      Cons2=Real1
      Cons3=Real0
      Go To 40
!

!** Mohr Coulomb_Chen
      5      Phira=Frict*Atan(1.0D0)/45.0D0
      Snphi=Sin(Phira)

      S3Ta=Sin(Three*Theta)
      Theta=Asin(S3Ta)/Three
      Thetat=(25.0/180.0)*Pi
      A1=1.432052062044D0
      A2=0.406941858374D0
      B1=0.544290524902D0
      B2=0.673903324498D0
      If (Abs(Theta).Le.Thetat) Then
         Kf=Cos(Theta)-Sin(Phira)*Sin(Theta)/Sqrt(Three)
      Else
         Sgn=Sign(One,Theta)
         Const1=A1+A2*Sgn*Sin(Phira)
         Const2=B1*Sgn+B2*Sin(Phira)
         Kf=Const1-Const2*S3Ta
      End If

      Cons1=Sin(Frict*Atan(1.0D0)/45.0D0)/Real3

      Alfa=1.0D0
      If (Abs(Theta).Le.Thetat) Then
         Dkf=-Sin(Theta)-Sin(Phira)*Cos(Theta)/Sqrt(Three)
      Else
         Dkf=-Three*Const2*Cos(Three*Theta)
      End If

      Cons2=Kf-Tan(Three*Theta)*Dkf
      Cons3=-Sqrt(Three)*Dkf/(Two*Cos(Three*Theta)*Sbar*Sbar)

      If (Abs(Theta).Le.Thetat) Then
         Cons2=Kf-Tan(Three*Theta)*Dkf
         Cons2=Alfa*Cons2
         Cons3=-Sqrt(Three)*Dkf/(Two*Cos(Three*Theta)*Sbar*Sbar)
         Cons3=Alfa*Cons3
      Else
         Cons2=Kf
         Cons2=Alfa*Cons2
         Cons3=Zero
      End If

      Go To 40


      40    Continue
      Do Istr1=1,Nstr1
         Avect(Istr1)=Cons1*Veca1(Istr1)+Cons2*Veca2(Istr1)+Cons3*Veca3(Istr1)
      End Do
      Return
      End

!********************************************************************************
!        *************** Get Fssi Mc Finished *************************
!********************************************************************************
